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(Received  20  August  1990;  accepted  23  October  1990) 

The  thermal  rate  coefficient  for  the  prototype  reaction  H  +  H3  —  H3  +  H  with  zero  total 
angular  momentum  is  calculated  by  summing,  averaging,  and  numerically  integrating  state-to- 
state  reaction  probabilities  calculated  by  time-independent  quantum-mechanical  scattering 
theory.  The  results  are  very  carefully  converged  with  respect  to  all  numerical  parameters  in 
order  to  provide  high-precision  benchmark  results  for  confirming  the  accuracy  of  new  methods 
and  testing  their  efficiency. 


I.  INTRODUCTION 

Quantum-mechanical  calculations  of  chemical  reaction 
rates  are  currently  practical  only  for  simple  systems,1-5  but 
they  can  be  used  to  test  approximate  methods  that  are  more 
generally  applicable.  As  new  methods  and  approaches  be¬ 
come  available,  it  is  useful  to  have  a  clearly  documented 
benchmark  that  can  be  used  for  testing  their  accuracy  and 
efficiency.  It  is  the  goal  of  this  and  the  following  paper*  to 
calculate  the  converged  quantum-dynamical  rate  coefficient 
for  one  well-defined  prototype  case  by  two  entirely  indepen¬ 
dent  methods  to  provide  such  a  benchmark.  The  case  chosen 
for  study  is  the  reaction  H  +  H,-H,  +  H,  where  the  atoms 
are  treated  as  distinguishable,  the  total  angular  momentum 
is  zero,  and  the  potential-energy  surface  is  given  by  a  double 
many-bodv  expansion  presented7  previously.  The  method 
used  in  this  paper  involves  time-dependent  scattering  theory 
to  compute  full  sets  of  state-to-state  transition  probabilities 
at  a  series  of  total  energies,  followed  by  summing  over  final 
states  corresponding  to  reaction,  and  averaging  over  initial 
states  and  energies  according  to  a  Boltzmann  distribution. 

Section  II  presents  the  theoretical  formulation.  Section 
III  describes  the  calculations  and  presents  the  benchmark 
results  and  convergence  studies.  Section  IV  gives  a  summa¬ 
rizing  conclusion. 

II.  THEORY 

Since  total  angular  momentum  ( J )  is  rigorously  con¬ 
served,  we  can  separately  evaluate  the  contribution  of  binary 
collisions  with  each  value  of  the  total  angular  momentum  to 
the  thermal  reaction  rate  coefficient.  The  J  =  0  contribution 
to  the  rate  coefficient  at  temperature  T may  be  written'' 

^  1  C  ' 

kJ  *  {T)  = - -  dEe  h:/kTN"(E),  (1) 

Jo 

where  fis  the  total  energy,  k  is  Boltzmann's  constant.  »t>*  is 
the  reactants’  partition  function  per  unit  volume  in  the  cen- 
ter-of-mass  frame,  and  XJ{  E)  is  the  cumulative  reaction 
probability'  (hereafter  referred  to  as  the  CRP).  The  CRP  is 
defined  as  the  sum  over  all  state-to-state  ( n  —  n ')  reactive 


transition  probabilities  PJana  „  ( E )  from  a  given  initial  chem¬ 
ical  arrangement  a  to  all  final  reactive  chemical  arrange¬ 
ments  a': 

lVJ(E)=  ^  P^an^E),  (2) 

o>  I 

n,n’ 

where  n  denotes  the  collection  of  all  initial  quantum 
numbers  (arrangement  a,  vibration  v,  rotation  j,  and  orbital 
/),  and  n'  denotes  the  set  of  final  ones.  (Note  that  a  and  a' 
are  also  specified  explicitly,  although  they  arc  contained  in  n 
and  n\  respectively. ) 

III.  CALCULATIONS  AND  CONVERGENCE  STUDIES 
A.  Calculations 

Converged  quantum-dynamics  calculations  were  car¬ 
ried  out  using  a  double  many-body  expansion7  of  the  poten¬ 
tial-energy  surface.  The  high  accuracy  of  this  energy  surface 
has  recently  been  verified  by  new  state-of-the-art  electronic 
structure  calculations."’  The  quantum-dynamics  calcula¬ 
tions  were  carried  out  by  the  generalized  Newton  variational 
principle  for  the  T matrix  involving  a  Lebesgue  square  inte¬ 
grate  ( l/'-)  expansion  of  the  reactive  amplitude  den¬ 
sity."  ’’  A  formulation  in  terms  of  the  T  matrix  with  com¬ 
plex  boundary  conditions  for  the  distorted-wave  radial 
functions  and  the  radial  Green's  functions  was  used.  Details 
of  the  basis  sets  and  numerical  procedures  for  this  kind  of 
calculation  are  given  elsewhere.'" 13  and  the  parameters 
used  for  the  present  study  are  based  on  previous  determina¬ 
tions14  of  converged  parameters  for  state-to-state  transition 
probabilities  for  this  reaction. 

We  carried  out  calculations  for  J  =  0  at  262  total  ener¬ 
gies  in  the  range  0.27-1. 66  eV.  Equation  ( 1 )  for  k  "  was  eval¬ 
uated  by  fitting  the  calculated  values  of  N"{E)  with  cubic 
splines  and  using  repeated  Gauss-Legendre  quadrature  to 
perform  the  integration.  Plots  of  e  '  h/kTN"(E)  for  200. 600, 
and  iOOO  K.  shown  in  Fig.  1,  exhibit  the  shape  typical  of  a 
thermal  distribution  of  reactive  collisions.  The  prominent 
shoulder  in  the  1000  K  plot  is  due  to  a  marked  step-like 
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FIG  I.  Bolfzmann-weighted cumufarive reaction  probability  forJ  =  0.  (a) 

:00K;  (b)  600  K;  (c)loQOK. 


feature  in  the  CRP  at  0.978  eV,  which  has  been  related  to  the 
threshold  of  a  linear  triatomic  quantized  transition  state.15 

The  partition  function  d>*  for  distinguishable-atom  H, 
without  nuclear  spin  is  given  by  the  expression 

(Inn'kl"'2 


<J)^  — 


V 


1(2/ 


+  l)e 


■  t../kT 


(3) 


where  n  is  the  reduced  mass,  h  is  Plank  s  constant,  and  e„  y  is 
the  vibrational-rotational  energy  of  the  reactant  H-,  diatom. 
The  conversion  factor  used  to  convert  from  atomic  units  to 
cm3  molecule  ~ 1  s'1  is  6.126  17.  The  rate  constants  for 
eight  temperatures  in  the  range  200-1000  K,  calculated  ac¬ 
cording  to  Eq.  ( 1 ),  are  reported  in  Table  I. 


B.  Convergence  studies 

A  careful  analysis  was  performed  to  determine  the  lim¬ 
its  of  accuracy  in  the  rate  constants  shown  in  Table  I.  There 
are  three  possible  sources  of  error  associated  win  the  calcu¬ 
lation:  (i)  convergence  of  the  quantum-dynamical  calcula¬ 
tions,  (ii)  numerical  integration  of  the  Boltzmann-weighted 
CRP,  and  (iii)  truncation  of  the  integral  in  Eq.  ( 1)  to  the 
energy  range  over  which  quantal  calculations  were  per¬ 
formed.  The  estimated  possible  error  associated  with  these 
sources  is  summarized  in  Table  II.  The  rest  of  this  section 
describes  how  we  obtained  these  error  estimates. 

The  first  source  of  error  is  from  convergence  of  the 
quantal  results.  This  was  evaluated  by  performing  calcula¬ 
tions  with  three  different  parameter  sets,  shown  in  Table  III. 
A  detailed  explanation  of  most  of  the  parameters  in  Table  III 
can  be  found  elsewhere.1213  Five  of  the  parameters  listed  in 
Table  III  have  not  been  discussed  fully  in  previous  publica¬ 
tions.  These  are  the  screening  parameters  IEPSTS, 
IEPSRD.  IEPSBS.  IEPSWM,  and  IEPSBM  which  simplify 
the  reactive  scattering  calculation  by  allowing  certain  quan¬ 
tities  to  be  set  to  precisely  zero  when  their  magnitude  is  very 
small.  Background  discussion  pertaining  to  the  explanations 
of  these  parameters  is  given  in  Ref.  1 2;  a  detailed  explanation 
follows  here. 

( 1 )  The  minimum  bond  distance  for  which  the  vibra¬ 
tional  wave  function  with  the  highest  value  of  v  and  j  —  0 
exceeds  10~IEPSTS,  and  the  maximum  bond  distance  for 
which  the  vibrational  wave  function  with  the  highest  value 
of  /  and  the  highest  v  for  that  j  exceeds  10  * IEPSTS  are  deter¬ 
mined.  In  calculating  interarrangement  (a=£a')  integrals  of 
the  coupling  potential,  angles  between  the  Jacobi  transla¬ 
tional  coordinates  R „  and  Rn  that  correspond  to  vibrational 
distances  outside  this  range  are  excluded  from  the  inner  an- 


TABLE  I.  Rale  constams  in  cm'  molecule  1  s 


T  (K) 

k  " 

:oo 

6.428  <.  10 

300 

8.505x10  “ 

400 

1.291x10  " 

MX) 

1.988x10  " 

700 

4.250X10  15 

1000 

1.578x10  14 
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TABLE  11.  Percent  error  in  k 

#  - 

vTiallKJIO. 

i  lumar.  ana  ocnwenne: 

i  nermai  reaction  rates  i 

Source  of  error 

200  K 

300  K 

400  K 

t>00K 

700  K, 

1000  K 

Convergence  of  quantal  CRPs 

0.028- 

0.024 

0.024 

0.022 

0.022 

0.021 

Point  density 

0.0059 

0.0027 

0.0029 

00025 

0.0029 

0.0034 

Truncation  of  integral 

0 

0 

0 

0 

0 

0.0071 

Total 

0.034 

0.027 

0.027 

0025 

0.025 

0.032 

-  This  should  be  interpreted  as  a  relative  error  of  2.8  x  10  4 ,  i.e.,  2.8  x  10 ' 1  %  ■ 


gular  quadrature  loop.  If  this  procedure  yields  a  negative 
value  for  the  minimum,  0  is  used  instead. 

(2)  In  performing  quadratures  at  small  values  of  the 
Jacobi  translational  coordinate,  only  distances  for  which  all 
radial  functions  (of  the  regular  solutions  to  the  distortion 
potential  problem  and  of  the  half-integrated  Green’s  func¬ 
tions1’)  have  magnitudes  greater  than  10“ lEPSRD  times  their 
maximum  value  are  included  in  the  integrals. 

( 3 )  When  Gaussians  are  used  for  the  translational  basis, 
as  is  the  case  here,  they  are  set  to  precisely  0  when  their  value 
falls  below  10  ' ,E,>SBS  times  their  maximum  value. 

(4)  Only  values  of  R„  and  Rn  +  ,  for  which  at  least  one 
matrix  element  of  W12  has  magnitude  greater  than 
10  " IEPSWM  are  included  in  the  exchange  integral. 

( 5 )  Only  values  of  R„  and  R„  +  1  for  which  at  least  one 
matrix  element  of  S'2  has  magnitude  greater  than 
10  " IEPSBM  are  included  in  the  exchange  integral. 

Set  A  in  Table  III  is  the  parameter  set  used  to  calculate 


TABLE  III.  Sets  of  parameter  values  for  GNVP  calculations. 


Parameter 

Explanation  Set  A 

Setfi 

SetC 

(u  —  0) 

a 

13 

14 

b 

jm.,  («=  1) 

a 

12 

13 

b 

Jmu,  <f  =  2) 

a 

n 

12 

b 

;m„.  U’=  3) 

a 

ii 

12 

b 

;«>  (u  =  4) 

a 

10 

11 

b 

U.,  («  =  5) 

a 

8 

m 

a 

10 

12 

b 

R  f  (a.u.) 

a 

2.047 

1.880 

b 

A  (a.u.) 

a 

0.335 

0.325 

b 

R‘;,  (a.u.) 

a 

5.062 

5.455 

b 

c 

a 

1.4 

1.3 

b 

N(HO) 

a 

75 

90 

b 

R„  i  (a.u.) 

a 

0.336 

0.325 

0.226 

R„  n.h  (Ml 

a 

12.909 

13.100 

15.490 

a 

10 

12 

b 

y 

a 

60 

70 

b 

a 

60 

70 

b 

V 

a 

12 

13 

b 

a 

14 

15 

b 

IEPSTS 

c 

9 

b 

20 

1EPSRD 

c 

20 

b 

30 

IEPSBS 

c 

20 

b 

30 

IEPSWM 

c 

20 

b 

30 

IEPSBM 

c 

20 

b 

30 

'See  Refs.  12  and  13  for  an  explanation  of  these  parameters. 
"Same  as  for  set  A. 

'  See  Sec.  Ill  B  for  an  explanation  of  these  parameters. 


the  CRPs.  Set  B  constitutes  an  increase  with  respect  to  set  A 
in  the  orders  of  quadratures,  the  size  of  the  basis  set,  the 
range  and  overlap  of  the  basis  functions,  and  the  range  and 
grid  point  density  of  the  finite  difference  grid.  In  set  C  we 
varied  the  parameters  not  varied  in  set  B,  and  two  of  the 
parameters  which  were  varied  only  minimal  amounts  in  set 
B  are  varied  more  significantly. 

The  difference  in  the  CRPs  calculated  with  parameter 
sets  A  and  B  is  shown  as  a  function  of  energy  in  Fig.  2.  For 
the  range  of  energies  which  contributes  significantly  to  the 
rate  constant  for  the  temperatures  studied  here  (£<  1 .2  eV), 
the  differences  in  the  values  of  N'\E)  calculated  with  the 
two  parameter  sets  are  very  small,  less  than  0.04%.  Error  in 
the  quantum-dynamical  calculations  was  parametrized  with 
the  straight  lines  shown  in  the  figure,  which  constitute  an 
estimate  of  the  error  in  the  CRP  calculated  with  parameter 
set  A.  In  order  to  represent  a  bound  on  the  error  in  k°  due  to 
lack  of  convergence  in  the  quantal  calculations,  a  value  for 
the  rate  constant  was  calculated  with  all  of  the  CRPs  in¬ 
creased  by  an  amount  corresponding  to  the  linear  parametri- 
zation  of  error  in  Fig.  2.  The  difference  between  the  original 
and  the  recalculated  rate  constants,  shown  in  Table  II.  is  less 
than  or  equal  to  0.028%  for  all  temperatures  studied. 

The  parameters  varied  in  set  C  have  less  effect  on  the 
results,  and  we  can  illustrate  this  by  considering  the  same  set 
of  energies  for  which  sets  A  and  B  were  compared.  The  lar¬ 
gest  error  in  the  CRP  calculated  with  set  C  relative  to  that 


FIG.  2.  Relative  percent  difference  in  the  CRP  as  calculated  with  the  Pr0- 
duction  parameter  set  A  and  the  convergence-check  parameter  set  6  The 
circles  are  data  points,  and  the  straight  lines  represent  the  parametrium'" 
used  to  estimate  the  possible  er~cr  in  A  "  due  to  error  in  the  quantal  CR 
caused  by  ihe  parameters  varied  in  set  B. 
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calculated  with  set  at  is  3  x  10  ~ 2  %  (at  0.4  eV),  and  there  is 
no  correlation  of  larger  errors  with  higher  energies  (the  er¬ 
ror  at  1.66  eV  is  only  5x  10  4%).  Thus  the  calculations 
performed  with  set  A  are  extremely  well  converged  with  re¬ 
spect  to  the  parameters  varied  in  set  C. 

The  second  source  of  error  is  that  associated  with  the 
numerical  integration  of  the  Boltzmann-weighted  CRP  be¬ 
cause  the  number  of  data  points  is  finite.  In  order  to  evaluate 
this  source  of  error,  we  calculated  the  effect  on  the  rate  con¬ 
stant  of  using  only  subsets  of  the  full  set  of  262  energies  in  the 
calculation.  To  do  this,  we  omitted  one  of  every  three  ener¬ 
gies  in  the  full  data  set.  Since  the  cubic  spline  fit  is  expected  to 
be  more  sensitive  to  the  omission  of  a  data  point  at  a  region  of 
hich  curvature  in  the  function  e  ~  E/kTN°(.E),  the  points  to 
be  omitted  were  chosen  in  three  different  ways  in  order  to 
provide  a  comprehensive  view  of  their  effect.  In  particular, 
e'  ery  third  point  was  omitted  beginning  with  either  the  first, 
t:.e  second,  or  the  third  point.  The  relative  difference 
between  the  rate  constant  calculated  with  points  omitted  and 
thit  calculated  with  the  fuli  data  set  is  less  than  or  equal  to 
0  )059%  for  all  temperatures  studied,  as  reported  in  Table  II 
( •  e  error  reported  for  each  temperature  is  the  largest  of  the 
three  differences  obtained  by  the  above  procedure). 

The  third  source  of  error  is  the  finite  upper  bound  of  the 
imegration.  This  was  studied  by  truncating  the  integration 
at  a  smaller  upper  bound  and  comparing  it  to  the  results 
detained  with  the  full  energy  range  of  0.27-1.66  eV  (since 
l  be  zero-point  energy  of  the  H:  vibrational  motion  is  0.27  eV, 
the  CRP  is  nonzero  only  at  energies  greater  than  0.27  eV, 
and  this  energy  can  be  used  as  the  lower  bound  for  the  inte¬ 
gration  without  causing  truncation  error) .  The  relative  error 
in  the  rate  constant  is  insensitive  (to  six  decimal  places)  to 
the  upper  bound  of  the  integration  until  the  temperature 
exceeds  700  K.  An  estimate  of  the  error  was  made  by  ap¬ 
proximating  iV°(£)  as  a  linear  function  at  energies  above 
1.66  eV  (see  Fig.  3)  based  on  an  appropriate  fit  to  the  enve¬ 
lope  of  the  oscillatory  IV'\E)  in  the  range  1.4-1.66  eV.  The 
possible  error  in  the  rate  constant  is  then  equated  to  the  inte¬ 
gral  of  the  product  of  the  linear  function  and  the  appropriate 
Boltzmann-weightmg  factor  from  1.66  eV  to  infinity.  The 


Energy  ieV) 


FIG.  3.  Cumulative  reaction  probability  lor  7  =  0.  The  straight  line  is  the 
linear  approximation  to  .V"(  E)  used  to  estimate  the  contribution  to  the  rate 
constant  from  energies  above  I  66  eV,  the  highest  energy  for  which  we  cal¬ 
culated  converged  quantal  results. 


error  thus  calculated  is  less  than  10  4  %  for  T<,7 00  K,  and 
is  7  X  10  '*%  at  1000  K. 

Finally,  an  upper  bound  on  the  total  possible  error  in  the 
rate  constant  was  estimated  by  summing  the  individual  error 
estimates.  As  can  be  seen  in  Table  II,  the  totai  error  for  all  of 
the  temperatures  studied  is  less  than  or  equal  to  0.034%. 
Thus  these  values  should  serve  as  definitive  benchmarks  for 
t  -.ting  new  and/or  approximate  procedures.  One  such  use, 
c  .'.ecking  the  stability  of  results  obtained  by  a  time-depen- 
ent  theoretical  method,  the  flux  autocorrelation  meth¬ 
od,6816  is  presented  in  the  following  article.  As  we  shall  see,6 
the  two  calculations  differ  by  at  most  0.022%  for  the  tem¬ 
peratures  common  to  both  studies.  This  value  is  well  within 
the  tolerances  for  the  calculations  reported  here  and  sug¬ 
gests  that  these  tolerances  are  generous.  The  excellent  agree¬ 
ment  confirms  the  accuracy  of  both  calculations,  since  the 
methods  used  are  completely  different. 

IV.  CONCLUSION 

The  c.,  ,tinguishable-atom  rate  constant  for  the  reaction 
H  +  H;~  il:  +  H,  where  the  total  angular  momentum  is 
zero,  was  llculated  by  a  method  that  is  exact  to  within  the 
limits  imposed  by  the  Bom-Oppenheimer  approximation 
and  the  best  available  potential-energy  surface,  the  double- 
many-body-expansion  potential-energy  function  of  Ref.  7. 
Careful  analysis  places  an  estimated  bound  on  the  error  in 
these  results  of  at  most  0.034%  for  temperatures  in  the  range 
200-1000  K.  Agreement  with  rate  constants  calculated6  by 
the  flux  autocorrelation  method  to  within  these  tolerances 
confirms  the  accuracy  of  both  calculations. 
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